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Zero Temperature Thermodynamics of Asymmetric Fermi Gases at Unitarity 
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The equation of state of a dilute two-component asymmetric Fermi gas at unitarity is subject to 
strong constraints, which affect the spatial density profiles in atomic traps. These constraints require 
the existence of at least one non-trivial partially polarized (asymmetric) phase. We determine the 
relation between the structure of the spatial density profiles and the T = equation of state, based 
on the most accurate theoretical predictions available. We also show how the equation of state can 
be determined from experimental observations. 

PACS numbers: 03.75.Ss 
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We consider the T = thermodynamics of a dilute 
asymmetric Fermi gas comprising two species of equal 
mass with s-wave interactions at unitarity. This has 
been recently realized in ^Li experiments [J 0, S 0] ■ We 
shall discuss the phase structure in the microcanonical 
and grand-canonical ensembles, and its manifestation in 
cold atomic traps using the local density approximation 
(LDA). The theoretical treatment of the grand-canonical 
ensemble is much simpler than the microcanonical ensem- 
ble as it consists of only pure phases. We discuss here the 
most general model-independent equation of state satis- 
fying known constraints. For model-dependent analyses 
see UMM- 

a. Phase structure: We show the main defining fea- 
tures of a grand-canonical phase diagram in Fig. [TJ The 
two species are labelled "a" and "5". The symmetry 
a ^ b allows us to consider only the region below the 
Ma = Mb where the locally averaged number densi- 
ties and chemical potentials satisfy n(, ^ and /if, ^ /Iq 
respectively. The asymmetry of the system may thus be 
characterized by the dimcnsionless ratios: 



nbha < 1, 



(1) 



Note that only x measures a physical asymmetry. There 
are four distinct regions: Vac — the vacuum, Nq — the 
fully polarized phases {x = 0) comprising only species 
"a", PPq — partially polarized phase(s) (0 < a; < 1), and 
SF — the fully paired symmetric superfluid phase (a; = l). 

We shall discuss only two phase transitions: one be- 
tween the fully polarized phase N^ (where a; = 0) and a 
partially polarized phase PPa (0 < n;, < tIq), and another 
between a (possibly different) partially polarized phase 
and the symmetric fully paired phase SF (where a: = 1). 
At unitarity, phase transitions occur along rays charac- 
terized solely by their slope yx- The two transitions we 
shall discuss are thus described by the two universal pa- 
rameters j/o and yi, which naturally satisfy j/o ^ yi. A 
major point of this paper is to place an upper bound Yq 
on yo (?/o ^ ill), lower bound Yi on yi (Yi ^ yi), and to 
show Y{) < Yi, which implies that the inequality ?/o < Vi 
is strict. This directly implies the existence and stability 
of one or more partially polarized phase(s) PP g. P ossible 
phases in the region PPa include LOFF states states 
with deformed Fermi surfaces [ll| , and p-wave superfluid 
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FIG. 1: Grand-canonical phase diagram of a two-component 
Fermi gas at unitarity and T = 0. Various phases are sepa- 
rated by phase transitions along the straight lines extending 
from the origin with constant slopes y^- The dotted line fol- 
lows the sequence of phases in a sample trap. 



states [l2|. If several of these states exists and are stable, 
the corresponding phase transitions will be characterized 
by additional universal parameters y^- 

h. Functional forms of thermodynamic potentials: 
At unitarity, the energy density £{na,nf,) and the pres- 
sure 'P{iJia, tib) have the following form: 



£{na,nb) = fa [nag{x)] 



5/3 



(67r2)2/3?i2 
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(2a) 



(2b) 



Note that g{x) = f^/^{x), where f{x) was introduced 
in 'Sij: The use oi g{x) rather than f{x) significantly sim- 
plifies the formalism [9]. The T = thermodynamic 
properties of the system are completely determined by 
the functional form of g{x) or h{y). The number densi- 
ties and chemical potentials are simply = dV /d^a,b 
and /ia,6 = dE/duafi respectively. As we show here, the 
functions h{y) and g{x) are tightly constrained by cur- 
rent Monte Carlo simulations, analytic calculations, and 
experimental results. The energy density and pressure 
are related via the Legendre transform 9] : 



Villa, Mb) = Ma^la + M6"-b ~ £{na, rib) 



^£{na,nb). (3) 



c. Physical constraints: The thermodynamics of 
three phases are known. The vacuum has vanishing pres- 
sure Vva.c — 0, the fuUy polarized phase Nq is a free Fermi 

gas, 

VfgM = iPt^'J^ (4) 

and the pressure of the fully paired phase SF is symmetric 
in fia and fib, and is described by a single parameter ^, 

^ / X 4 /5 5/2 fia±fib .^s 

VsF{ti+) = 5^M+ where /z± = . (5) 

These provide the limiting forms for h{y) and limiting 
values of g[x), [see Eqs. (|6a| and ([7a| below]. 

The phase transition at y = tjq defines the border of 
the region where fib = J/oMa is tuned to keep species "6" 
out of the system. Since the interspecies interaction is 
attractive, the critical fib must be negative yo < 0. We 
will provide an upper bound Yq (yo ^Yq). 

Note that 'Psf{iJ'+) depends only on the average chem- 
ical potential This inscnsitivity to the chemical po- 
tential difference fj,- is due to the existence of an energy 
gap A in the spectrum. The phase transition at y = yi 
marks the line where fi- becomes large enough to break 
the superfluid pairs. In the phase SF, /x_ is constrained 
by the size of the physical gap fi- ^ A l8| . This provides 
a lower bound Yi ^ yi, see below and [9| . 

(Recall from Eq. U]) that we are only considering re- 
gions where < ^_ = Im-D- 

If no stable partially polarized phase exists, then the 
region PFq will vanish, being compressed into a single 
first-order transition line where pressure equilibrium is 
established VpcitJ-a) = 'Psf(m+) [1, Hi- This would occur 

y = Vc = (2^)'^/^ — 1 = yo = yi, and would imply that 
Yq = Yi. We argue below that Yq is strictly less than Yi , 
and therefore rule out this possibility. 

Finally, thermodynamic stability requires that the 
pressure and energy density are convex functions, which 
implies that g{x) and h{y) are also convex Q. The con- 
straints on h{y) are 



h"[y) ^ 0, and (6b) 

yo^Yo<yc<Yi^yi^ 1. (6c) 

The corresponding constraints on g{x) are 

5(0) = 1, 5(1) = (20'/', (7a) 

g"{x) ^ 0, and (7b) 



g'{0)^Yo, g'{l)e[g{l)(l + Y-^)-\g(l)/2]. (7c) 

Equation ((6a|) follows directly from Eqs. \2h\ . ([4]), and 
(O, Eq. ([7a|) follows from Eq. Pa|) . and the interval 
in Eq. ([7c)l follows from the properties of the Legendre 
transform and Eq. (pc)) 0. 
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FIG. 2: Example of a function h(y) and the corresponding 
function g{x) shown as thick lines. Maxwell's construction 
for phase coexistence leads to a linear g{x) for x £ (0.5, 1.0), 
interpolating between the two pure phases shown with lighter 
lines. This corresponds to the kink (first-order phase transi- 
tion) at 1/ = 2/1 in h{y). Various other sample functions are 
lightly sketched within the allowed (dotted) triangular region. 

d. Parameters: For Fig. [2] we used 

C 0.42(1), Yq w -0.5, Yi = -0.09(3), yi = 0.05. 

We obtain estimates for Yi and ^ from Monte Carlo 
data [13, [3, The latest Monte Carlo estimates 

for the symmetric systems give ^ — 0.42(1) [ij, fl5| 
and A/ep = 0.504(24) [l^, where ep is the Fermi en- 
ergy of the free gas with the same density. This gives 
j/c K, —0.099(15), and the constraint fi- ^_^A gives 
n = (e - A/£p)/(,e + A/ep) = -0.09(3) g Since 
within the statistical errors Yi ~ ?/c, the possibility of an 
empty PPq region at unitarity cannot yet be ruled out 
by this value of Yi , as was noted earlier by Cohen Q . 

We now estimate Yq- Let eo be the energy required 
to add one particle 6 to a fully polarized gas of density 
Ua- Consider adding a large, but infinitesimal amount 
of 6, 1 ^ iVb <C Na- In the thermodynamic limit, the 
required energy per particle will be the critical chemical 

potential fib — ar?J^ g' {0) defining the transition j/o- If 
the added particles repel, the energy will be NbC^, and 
fib = eo- If they bind, the additional binding energy must 
be included, giving fib < cq. In this way, e^) provides a 
bound for fib and 5'(0) = yo ^Yq — eo/(ana^ ). 

Consider adding a single b fermion, with coordinate Tq, 
to a system of Na a fermions with coordinates r„. Let 
Tnm = |rn — | . The wavc function for the b fermion in 
the background of fixed a sources is 

(/) ro; {r„} = > , 8) 
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and it satisfies the zero-range interaction boundary con- 
ditions if the following Na conditions are met: 



1 



= 0. 



(9) 



For uniform distributions where the lowest state has con- 
stant A„i = A, approximating the sum as an integral 
gives K—a^^ — Aima/K^. This continuum approximation 
is not very accurate in the unitary limit, since k is compa- 
rable to the inverse interparticle spacing and /ua ~ 47r. 
To estimate corrections, the equations can be solved for 
various lattice configurations. We find that k deviates 
from the continuum result by a factor of 0.84(3) for sim- 
ple lattice configurations and perturbations (see Q for 
details). We now estimate the energy of the system us- 
ing the wave function 



^{ro;{rn}) = $sD({r„})(/>(ro; {r„}). 



(10) 



where $sd is a Slater determinant for a free Fermi gas 
and obtain eo ~ —h^K^/m Q: 



■iTrfi^naa/m, 
eo ~ -0.71(5)ft2(47rna)^/^/TO, 
{ma?), 



if a 
if a 
if a 



±oo, 
0+. 



(11) 



Note that this result interpolates between the correct 
leading order BEC and BCS results. This estimate as- 
sumes that the fluctuations of the number density na(r) 
on a scale of the order 1/k affect k very little. The result 
is consistent with this assumption, as discussed in 
The constraint at unitarity is thus 



lo « -0.54(4) <Vc = -0.099(15). 



(12) 



If Yq is strictly less than t/c, then convexity in g{x) and 
h{y) implies yc < Yi (see Fig. [5]). 

e. Trap profiles: For large systems with a slowly 
varying confining potential, gradient terms may be ne- 
glected, and the LDA employed to determine the den- 
sity distribution by introducing spatially varying effective 
chemical potentials: 



(13) 



Lagrange multipliers Ao,& fix the total particle numbers 
Na and Ni,. The LDA may be inaccurate near phase 
boundaries where the densities change rapidly. The gra- 
dient terms will smear out these transition regions and 
provide an additional surface tension proportional to the 
local curvature p^ . 

In the LDA, the density profile may be constructed 
from the local densities ria^b using Eq. (|13p (explicit for- 
mulae are provided in Q). The dotted line in Fig. [1] 
shows the sequence of phases contained in a sample trap. 
Since 2fi^ — Xa — is fixed, traps contain the se- 
quence of phases encountered along a 45° line through 
such a diagram. In this example, the center of the cloud 
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FIG. 3: Measured transition radii from Ref. The upper 
plot shows the normalized data (crosses) on top of data gen- 
erated from several randomly generated functions h{y). The 
lower plot shows the parameter 7 defined in Rel. (|14[) . 



(F(0) — 0) is in the SF phase at the cross. The phase 
transitions will occur for y = yi, y = yo, and y = —00 
for l/(Ri)==Fi, F(Ro) = Vb, and y(R^^e) = Kac respec- 
tively, with Rq < Ri < i?vac- As noted above, additional 
phase transitions may exist between Rq and 

/. Experiments: Accurate measurements of the den- 
sity profiles would allow for a complete extraction of 
h{y) and g{x). For example, using x = nb(R)/na(R) 
and the expressions for fj,a,b, we have g'^/^{x)g'{x) = 

[Xb — V{'R)]/[ana^^ (R)] from which g{x) may be ex- 
tracted using the boundary condition g{0) = 1 |9|. 

For harmonic traps, the locations of the main phase 
R 

vac 



transitions, i?vac oc VKac, Rq oc Vn), and i?i cx vFT, are 
completely determined by the Lagrange multipliers Aa^^, 
and the universal numbers yo and yi. Within the LDA, 
we obtain the following model independent relationship, 
to be compared with the recent MIT data 3 (see Fig.[^: 



7 = 



1-2/1 
1-2/0 



^vac ^0 
^vac -"-1 



0.70(5). 



(14) 



To extract more information, one must consider a specific 
functional form for h(y) and g(x). We have analyzed a 
large sample of allowed functions h(y) and g{x), a few of 
which are sketched in Figs. [2] and [31 We find that the 
total polarization P = {Na — Nb)/{Na + Nb) is quite in- 
sensitive to the functional form. However, the critical po- 
larization Pc — where the innermost phase transition ap- 
proaches the center of the trap i?i = — is quite sensitive 
to yi. The MIT experiments measure P^ = 0.70(5). 
If yi = 0, then one obtains Pc > 0.80 . . .0.85. However, 
if one considers yi « 0.05 {g'{l) ~ 0.04), then values 
of Pc ~ 0.7 and smaller emerge, compatible with those 
measured in [ll,i]. Using Eq. nil), this gives a value of 
yo = g'{0) ~ -0.4. Our estimate for Yq w -0.54(4) is 
consistent with this extracted experimental value within 
existing uncertainties. 



4 



Within the Eagles-Leggett extension of the BCS 
model [13], one obtains the values yo — 0, yc — 0.105, 
and yi = 0.107 (see 0, H), which would correspond to 
a parameter 7 = 0.893, as opposed to the 7 = 0.70(5) 
extracted from experiment. At the same time, the spa- 
tial layer for the PPa region would be very thin, namely 
(i?g - Rl)/{Rl^^ -Rl) = l-j = 0.107, compared with 
our estimate 1 — 7 « 0.30(5). 

Our analysis is strictly valid only at T = 0. The devi- 
ations in Fig. [3] are most likely finite temperature effects. 
The regions of the phase diagram most sensitive to T ^ 
are those with small fi, thus, the transition radii in traps 
with small asymmetry will be most affected. For large 
polarizations, the temperature should not affect the SF 
phase which is gapped, but will affect phases with zero 
or small gap excitations. This could alter the values of 
yi slightly and yo significantly. We thus caution against 
taking the extracted numbers in this section too seriously 
until a similar finite temperature analysis is presented. 

Since Yi < (as determined from the value of the pair- 
ing gap), a positive t/i suggests a first-order phase tran- 
sition out of the SF phase for a critical < A. This 
conclusion is also consistent with the form of the quasi- 
particle spectrum computed in [isj . which has a mini- 
mum at a finite momentum. 

In conclusion, we have shown that thermodynamic con- 



straints, accurate Monte Carlo simulations, analytic esti- 
mates, and experimental data place tight constraints on 
the equation of state of the asymmetric T = unitary 
gas. These constraints imply that there exists a region 
where one or more nontrivial partially polarized phases 
exist. These phases likely exhibit very interesting mi- 
croscopic physics. In particular, any ungapped polarized 
phase is unstable towards the formation of a state with 
two symbiotic superfluids at T = [ll|. The tight con- 
straints on the forms of g{x) and h{y) we present will 
help locate these novel phases. 

The authors thank A. Schwcnk and M. W. Zwierlein 
for comments, D. T. Son for many useful discussions, 
M. W. Zwierlein et al. for providing the data in Fig. [31 
and the US Department of Energy for support under 
Grant No. DE-FG02-97ER41014. 

Notes added: Chevy [l8l | independently arrived at sim- 
ilar conclusions. The latest MIT analysis 19] agrees with 
our conclusion that the SF phase occupies the center of 
the trap and shows that the LDA is applicable. In a 
recent variational Monte Carlo study, Lobo et al. [l^] 
agree with our lower bound, obtaining Iq = —0.58(1) < 
yc = —0.099(15), and concluded that the transition at yi_ 
is first-order. This is consistent with our results: their 
fimction f{x) is very similar to our g^/'^{x) (see ^). 
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APPENDIX A: ADDITIONAL DETAILS 

In this appendix, we present some additional calcula- 
tional details. We discuss the properties of the Legendre 
transformation, thermodynamic stability and convexity 
of thermodynamic potentials, and the Maxwell construc- 
tion for mixed phases. We also provide details of the cal- 
culations whose results are presented in the main body. 



1. Properties of the Legendre Transformation 

Many of the thermodynamic constraints follow directly 
from properties of the Legendre transformation ([3]) . We 
explain some of these properties here, because there has 
been some confusion in the literature. 

The Legendre transform relates tangents in one en- 
semble to coordinates in the other. For example, a linear 
segment of g{x) over a finite interval in x will be mapped 
into a single point ?/, where h{y) has a kink. Straight seg- 
ments of g{x) arise from the Maxwell construction and 
indicate a phase coexistence (mixed phase) . For this rea- 
son, the grand-canonical phase diagram is much simpler 
than other ensembles as it consists solely of pure phases. 

First we present some relations. Starting with the 
definitions of the thermodynamic potential densities 
£{na,nb) and —'P{^,a, p^b) defined in we differentiate 
to find the densities 

"a = ^ = PMy)f'^[h{v) ~ yh:{y)l (Ala) 

nb = j^^ PMy)f/^h'[y), (Alb) 
and the chemical potentials 

fi^^ — = a[na9ix)f^^[g{x) - xg'ix)], (A2a) 

= ^ = ^K9{x)?^'9'{x). (A2b) 

From these relations and the definitions of the asymme- 
try parameters x = Ub/ua and y = ^bl y^a, we obtain the 



V = 



g{x) ~xg'{x)' 

h'jy) 
h{v) -yh'ijjY 



g{x) -xg'{x)' 



h{y) - yh'{y) 



(A3a) 



(A3b) 



The important geometric property of the Legendre trans- 
formation is that it maps tangents in one space to 
points in the other. Consider the grand-canonical en- 
semble described by the function V{^J,a, fJ-b)-^ The tan- 
gents {dfj_^V, df^i^V) — {ria, rib) to this function at a point 
ifJ-a, ^J■b) in the grand-canonical ensemble map directly 
to the coordinate {na,nb) in the microcanonical ensem- 
ble. Note that the Legendre transformation is symmet- 
ric: tangents in the microcanonical ensemble map back 
to points in the grand-canonical ensemble. 



2. Thermodynamic Stability: The Second Law 

The thermodynamic potentials follow from a strict 
minimization procedure over all possible states. If this is 
properly carried out, the potentials will be convex func- 
tions of their arguments. This convexity is the geometric 
manifestation of the second law. Locally, the requirement 
is that the Hessian of the potentials — the matrix of sec- 
ond partial derivatives — be positive semi-definite. One 
can easily show that this requirement also implies that 
both density and concentration sound modes are stable 
with corresponding real sound velocities. We shall eval- 
uate now the Hessian's for the thermodynamic potential 
densities in the case of a two species at unitarity and 
T = 0. We thus find that 



HiVh) cx 
HiSg) cx 



3(/i - yh')2 + 2y2hh" 3{hh' - yh'^) - 2yhh" 
3{hh' - yh'^) - 2yhh" 3h'^ + 2hh" 

4(5 - xg')2 + 6x2gg" ^{gg' - xg'^) - 6xgg"' 
4(5.g' - xg'^) - 6xgg" ig'^ + 6gg" 

10/ - Uxf + 9x^f" 6/' - 9xf"' 
6/' - 9xf" 9/" 



The corresponding determinants are: 

det[H{rh)] cx h^h" 
det[Hi£g)] cx g^g" 
det[Hi£f)] cx 5//" - 2/'2. 



(A4) 
(A5) 
(A6) 



From these it is easy to see that the Hessians of V and 
£ are positive definite if h" ^ and /i ^ 0, or g" ^ 
and 5^0. If one were to use the parametrization of the 



^ Strictly speaking the grand-canonical potential is —V'P{fia,fJ'b), 
notice the minus sign, where V is the volume. 
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energy density E through / as suggested by Cohen P 
the situation is more compHcated, as one would have to 
satisfy the nonhnear differential constraint 

5//" - 2/'2 ^ 0, (A7) 

c/. Eq. (6) in Ref. 01 . 

Though not related by a strict Legendre transforma- 
tion, equations (jA3p show that the universal functions 
g{x) and h{y) are similarly related, hence Fig. [2] exhibits 
the same Maxwell construction properties. Note that the 
linear Maxwell construction only works with the func- 
tional form of g{x) as introduced in (Pa|) : if one uses 
the function f{x) — g^/^{x) introduced in [s], then the 
equivalent construction will involve nontrivial forms of 
f {x) that saturate the inequality (IA7P Q: 

f{x) = {A + Bxf/^. (A8) 

For this reason, the function g{x) is substantially simpler 
to discuss than f{x). By further requiring saturation of 
inequality (|A7|) throughout the entire interval, Cohen fsl 
obtained 

f{x)^{l + mf'^-l]xY'\ (A9) 

which describes the case of a vanishing partially polarized 
region PPa- 

Finally, we discuss the physical significance of kinks 
in £{n) which translate to flat regions of 7^(/i) in the 
grand-canonical ensemble. A kink in £{n) means that, 
for a range of tangents (chemical potentials), the energy 
of the ground state does not change. In other words, 
the system is insensitive to a range of chemical poten- 
tials. For example, the superfluid phase SF is insensitive 
to a range of chemical potential splitting /i_ because of 
the physical gap in the spectrum. The flat regions in 
the grand canonical ensemble represents the same physi- 
cal state (same densities) that are stable over a range of 
chemical potentials due to this gap. 

As discussed in the main text, the system may or may 
not respond to chemical potential differences strictly less 
than the gap /i_ < A, but will definitely respond when 
/i_ > A. This gave us the bound Yi on the transition 
parameter yi (see (jAlSp below). 

3. Maxwell Construction 

As an example, let us consider the Maxwell construc- 
tion for phase coexistence in the microcanonical ensemble 
with a single species "a" . The Maxwell construction for 
the case of two species at unitarity is somewhat more 
complicated, but, as we discussed above, involves the 
same type of linear construction if the function g{x) is 
used. This ensemble is constructed by minimizing the 
energy density £{n) over all phases p with fixed total 
particle number N = nV , and volume V: 

£{n) =m:m£p{N/V). (AlO) 



In the microcanonical ensemble, this procedure is slightly 
complicated by the possibility of phase coexistence. Con- 
sider two distinct pure phases p G {1,2} with given en- 
ergy densities £i{n) and £2{n)'. an example of two such 
phases is shown in Fig. [4] Minimizing (|A10[) over these 
phases separately produces an energy density £{n) that 
is not convex. The Maxwell construction proceeds by 
constructing a series of systems with fixed N by combin- 
ing a physical fraction z of the system in phase pi with 
density ni and the remaining fraction 1 — z in phase p2 
with density n2 subject to the constraint of fixed average 
density 

TV 

rt = — = ZTii + (1 — z)n2. (All) 

This leaves two degrees of freedom and defines a series of 
two-component mixed phases that must also be consid- 
ered in (|A10[) with energy density: 

£{n) = z£p,{ni) + (1 - z)£p^{n2). (A12) 

Minimizing £{n) — pn, where p is a, Lagrange multiplier, 
leads to the linear segment shown in bold in Fig. 21 
which is the convex hull of the energy densities. More 
specifically, by varying 71,1^2 and (and assuming that 
all derivatives exists), one obtains the slope 

_ gpi(rj-i) - £p^{n2) _ d£p^{ni) _ d£p^ (n2) 
ni — 712 dnx dn2 

(A13) 

This is simultaneously the slope of the line segment defin- 
ing the mixed phase, and the tangents to the two energy 
densities. In the case of the shown in Fig. [H the 
function has a cusp at the point of transition. In this 
case, d£p^{n2) / dn2 should be interpreted as the appro- 
priate value in the interval defined by the corresponding 
derivatives computed to the left and to the right of the 
kink. In (jA13p dSp^l dn2 may assume any of these values 
to satisfy the equilibrium condition. This construction 




n 



FIG. 4: Maxwell construction for two competing phases in the 
microcanonical ensemble p G {1,2} with given energy densi- 
ties £\ (n) and £2 (n) . The right phase has a gap as signified by 
the kink. The left phase has no gap. A first-order phase tran- 
sition connects the two phases. Minimizing HA10|) produces 
the thin solid curve, but this is not convex. The Maxwell 
construction amounts to finding the convex hull — the thick 
solid line — which restores the convexity of £{n) required by 
the second law of thermodynamics. 
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guarantees that the potential £{n) will be convex as re- 
quired by the second law of thermodynamics. 

The conditions of minimization are equivalent to the 
conditions that the chemical potential and pressure of the 
two coexisting phases be equal. Note that phase coexis- 
tence occurs only along the linear segment. Along this en- 
tire segment, the tangent is the same. Thus, the entire re- 
gion of phase coexistence is described by a single point /^c 
representing the phase transition in the grand-canonical 
ensemble. Furthermore, since the density changes dis- 
continuously from one side of the phase transition to the 
other, the pressure — though continuous — will have a kink 
at jjic, consistent with a first-order transition. 

For this reason, phase structures are much simpler 
when considered in the grand-canonical ensemble. In this 
ensemble, the phase diagram consists of only pure phases: 
all phase coexistence (mixtures) occur along first-order 
phase transitions. 



These constraints transform directly into constraints on 
the derivatives of g{x) through the dictionary (|A3p . For 
example, when a; = we have 



yo < >o 



^-g'(0)<yo, (A19) 



where we have used the fact that 17(0) = 1. Likewise at 
a; = 1 we have 



yi>Yi 



which gives a lower bound on g' (V). The upper bound on 
g' {V) is simply the condition that yi < 1 which we may 
impose by symmetry: 



yi < 1 



4. Constraints 

Here we present a few more details about the extrac- 
tion of the various constraints in the main text. We start 
with the bound Yq. Recall that the transition is defined 
by chemical potential /if, required to keep out species "6" . 
This is bounded by the energy gained by adding the sin- 
gle particle eo which we estimate: fit, < cq. We must re- 
late this to the chemical potential through the density 



of the free Fermi sea Ua = PfJ^a 



3/2 



an 



2/3 



Ma = K//3)2/3 



We now express the constraint < eo in terms 



of yo: 



, eo 

= — - —y^ 

H-a ana 



= Yo. 



(A14) 



The bound Yi follows from the constraint /i_ < A. To 
express this in terms of the data, we need to express the 
chemical potentials in terms of the normalization Fermi 
energy ep = /^+^ of a free gas of the same density n+ — 
na + nfc = 2(2m/iFG)3/2/(g^2). 



2m 



2/3 



(A15) 



The relationship between the density of the symmetric 
phase SF and the chemical potential /x+ is determined 
by the parameter ^ from 



n+ 



dVsFO^+) _ 213 3/ 



^3/2 



(A16) 



Thus, we have simply ep — tJ-+ /£,■ We now express the 
constraint: 



Solving for yi gives the constraint 



yi>Yi 



(A17) 



(A18) 



5. Density Profiles 

We consider the following functional form for g{x): 

g{x) = \ ^ ^ where x € [0, xt], ^^22) 
I c + where x S [xt, 1] 

The form of the function has been arbitrarily chosen be- 
tween and Xt (for simplicity, we chose a quadratic 
polynomial) with the correct intercept g(0) — 1 and 
slope (?'(0) — 2/0- From this curve, we proceed with 
the Maxwell construction by finding the line that passes 
through {x,g{x)) = (1,(7(1)) and that is tangent to the 
given curve g{xT) at the point xt- The resulting g{x) is 
shown with solid thick line in Fig. [2] and Fig. [T] 

Once this curve is established, one can extract hiif) 
using Eqs. (|A3p . The only complication arises when there 
are kinks or linear segments. In this case, there is a linear 
segment of g(x) for x £ [xt, !]• From (|A3p we see that, 
throughout this region, y and hijj) take on only a single 
value: this is the first-order phase transition . This will 
result in a kink in the function h{y) at this location. 

To compute the density profiles in a trap, we first 
parametrize the trap so as to establish the spatial varia- 
tions of the chemical potential. Let us consider a spher- 
ical harmonic trap V{R) oc and consider the coordi- 
nate R — R/Rvac where i?vac is the radius of the cloud. 
The effective chemical potentials are established by 
once the Lagrange multipliers Xa^b are chosen. Once this 
is done, the functional form of h{y) can be used to di- 
rectly map the position R to the densities using Eq. (|Aip . 
In Fig. [S] we show several trap density profiles with the 
sample function (|A22p . 



6. Extracting g(x) from Experiment 

In principle, once the radial density profiles have be 
measured to sufficient accuracy, one can extract the func- 
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0.25 




^ and A (see for example f^): 



FIG. 5: Density profiles na,b of the two species in a spherical 
harmonic trap as a function radius R = R/Rvac in units of the 
cloud radius i?vac for thermodynamic function (|A22|) . Density 
profiles are plotted for a fixed A+ — {Xa + Xt) /2 and a variety 
of chemical potential differences jj,- — X- ranging from A- = 
(fully paired ria ~ rib throughout the trap) to A_ = A+ (no 
superfluid core). The red solid lines are the majority species 
Ha while the black dotted lines are the minority species. The 
critical radii for the intermediate profile have been denoted 

^0,1- 



a/F = 0.5906, 



= 0.6864. 



(A24) 



These determine the properties of the superfluid phase 
SF. If we consider only homogeneous and isotropic 
phases, then there are three distinct phases: SF, Na b, 
and Nq.^ The SF phase is the usual symmetric BCS- 
BEC crossover phase, the PPa phase is a partially po- 
larized two-component Fermi liquid, and the Nq phase 
is a fully polarized single-component Fermi liquid. One 
of the shortcomings of the mean-field crossover models is 
that they neglect the Hartree-Fock contributions. 

These terms enter the energy density as Ehf ~ grianb 
where g is the coupling constant. In the limit of a short- 
range interaction, one takes the range of the interaction 
To — !■ to zero while holding the vacuum inverse s-wave 
scattering length fixed. This requires that the in- 
teraction strength be taken to zero g ~ Trro/m —^ 0. The 
densities na,b contain no singularities and so the Hartree- 
Fock contributions vanish. These contributions may be 
included in weak coupling by resumming particle-particle 
ladders to obtain £hf ~ 4:Tr fi^ an arib/m, but this proce- 
dure may not be extrapolated to unitarity. 

Thus, in mean-field, the partially polarized state sim- 
ply has the pressure of two independent free Fermi gases: 



^FG.(A'a,/ib) = |/3(Ma/'+/4^') 



(A25) 



Pressure equilibrium thus determines the phase transi- 
tions in this model and we have the two conditions 



^FG2(Ma,yf^^Mb) = ^Sf(M+) 



(A26a) 
(A26b) 



tional forms for g{x) and h{y). For example, an exper- 
iment can extract the density profiles na^b(R) (see for 
example [l^) for a known trapping potential V^(R). Us- 
ing relations ([13]) and (jA2bp . we have: 



g'/MR)]gW)] = 



Xb - VCR) 
aua (R) 



(A23) 



This first-order differential equation may be integrated 
with the boundary condition 17(0) = 1 to find the function 
g{x). Finally, Xb will have to be fit to the trap profile. 
Note, there are many other ways to extract the same 
information: we have simply chosen a simple method. 
Other methods may be less sensitive to experimental er- 
rors for example. We leave it to future work to perform 
this extraction and the accompanying error analysis. 



These imply trivially that j/g^^ = 0. Since the Hartree- 
Fock terms vanish, there is no interaction energy and 
the phase transition between N^ and PP^ occurs for 
/ib = 0. The SF/PPq. transition at yi is governed by 
condition (jA26bp . Using Eq. ([5]) and Hb — yoUa, we have 



9t-3/2 



1 



Hi 



MF\5/2 



(A27) 



This may be solved numerically using Eq. (jA24p to ob- 
tain yi^^ — 0.1067. To summarize: In the Eagles- 
Leggett mean-field crossover model with homogeneous 
and isotropic phases, we have: 



= ro^-^^ = o. 



r*^^ = -0.0750, 



^ = 0.1051, 
2/f ^ = 0.1067. 



(A28a) 
(A28b) 



7. Mean Field (Eagles-Leggett) Results 

In this section we consider the Eagles-Leggett mean- 
field model [l3|- In this model, one can easily calculate 



^ There is the possibility that other phases, like LOFF, compete 
with the partially polarized phase In the mean-field model 

it is unlikely that these phases will drastically alter the locations 
of the phase transitions. The LOFF-like regions are always very 
thin, and we neglect these possibilities here. 

^ The = limit at unitarity does not represent a singularity. 



9 



Obviously the lower bound is useless for our pur- 

poses. 



which arises from applying Hq to both the Slater deter- 
minant and (/'(ro; {fn})- This gives a vanishing contribu- 
tion, since the integral over momenta 



8. Calculation of eo 



The Hamiltonian describing a system of Na species 
"a" fermions interacting with one additional species "&" 
fermion is: 



H = H„ + Ho = ^ t„ + to + ^ V„,o , (A29) 



where Tfe are the corresponding kinetic energy operators 
and V„^o is the potential between the fermion of species 
"a" with coordinate r„ with the fermion of species "b" 
with coordinate tq. Recall that we use the variational 
wave function (jTU]) 



*(ro; {r„}) = $sD({r„})0(ro; {r„}), 



where 



(/)(ro;{r„}) = 



exp i-Kron) 



ron 



(A30) 



(A31) 



is an eigenfunction of the operator Hq 



Ho^(ro;{r„}) 



2m 



0(ro;{r„}). (A32) 



The Slater determinant is of course an eigenfunction 
of H,: 



Ha$SD({r„}) = EFG^SB{{rn})- 



(A33) 



Two additional terms arise from Hq acting on the prod- 
uct wave function. The first term is a cross-derivative 




0.0 0.1 0.2 0.3 0.4 

Maximum displacement (units of lattice spacing) 



m 



dk 

(2^ 



Na 



Na 



1=1 



(A34) 

is identically zero because of the symmetry of the Fermi 
sea. The second kind of term is due to Hq acting on 
(/)(ro; {fn}) and which pulls down derivatives of K({r„}) 
and A„({r„}): 



Ha(/)(ro;{r„}) = - 



0(i"o ;{!■„}) 



2m 

■ derivative corrections. 



(A35) 



In what follows, we neglect the derivative corrections, 
but, as we discuss below, we expect these to be small. 
Within this approximation, the expectation value of this 
Hamiltonian is 



(*|H|*) =EFG + ea=EFG 



2^2 



(A36) 



where Epo is the ground state energy of the system of 
Na fermions alone. This contribution arises by applying 
'^n=i -^n ^'^ Slater determinant alone. It is implied 
here that 4>{ro; {r„}) is normalized. 



We have checked explicitly that at unitarity for several 
simple lattice configurations (fee, bcc) and for configura- 
tions in which the positions of all the fermion species "a" 
were randomly changed from the ideal lattice configura- 
tions within the Fermi hole, the values of k varied by a 
few percent at most, see Fig. This analysis thus 

confirms our assumption that the fluctuations in a free 
Fermi gas of the number density na(r) do not affect in a 
noticeable manner the value of k. 



FIG. 6: Relative change in k for randomly perturbed bcc 
lattices. The abscissa shows the maximum deviations each 
lattice site was displaced from equilibrium as a fraction of the 
lattice spacing. Curves are shown for cubic lattices from 10 to 
14 sites per side. Each curve shows the maximum deviation 
over a sample of 30 random lattice configurations. 



For comparison, we replot the function q(x) from Fig. [2] 
and include the Monte Carlo data from [20] . The authors 
of Ref. 20j simulated a two-component polarized normal 
Fermi gas, which provides a variational bound on the en- 
ergy. We thank S. Giorgini for sending us their numerical 
results (2Q] . 
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FIG. 7: Monte Carlo variational upper bound on g{x) 
from |20l | plotted on top of the function g{x) from our Fig- 
ure[2] Note the agreement from small polarizations indicating 
that our estimate for Yq is consistent with their proper varia- 
tional bound. For larger polarizations, the true curve will lie 
below the results of Ref. [20| for two reasons: 1) The Maxwell 
construction for g{x) (see Fig. (2) of [20l |') and 2) The authors 
of Ref. [2^ considered only normal Fermi partially polarized 
states. As shown in fl^l, at T = 0, partially polarized states 
will be superfluid. This could noticeably lower the energy for 
substantial polarizations at unitarity. 



